#import brian_no_units
from brian import *
from brian.library.random_processes import *
from brian.library.synapses import  *


mS = msiemens

# Conductances (mS/cm^2)
GNa = 120 * mS/cm**2
GK_dr = 100 * mS/cm**2
GCa_NS = 14 * mS/cm**2
GCa_ND = .03 * mS/cm**2
GK_CaS = 3.136 * mS/cm**2 # canonical: 5
GK_CaD = 0.69 * mS/cm**2  # canonical: 1.1
GCa_L = 0.33 * mS/cm**2
gleak = 0.51 * mS/cm**2
GNapD = 0.1 *mS/cm**2

# Static parameters
C = 1 * uF/cm**2
gc = 0.1  * mS/cm**2      # coupling conductance (mS/cm^2)
p = 0.1
Kd = 0.2 * 10**-2*mole*dmetre**-3     # uM
f = 0.01      # percent free to bound Ca
alpha = 0.009 *mole/(amp*second)*um**-1# mol/C/um
kca = 2 * ms**-1      # Ca removal rate



# Half Activation Vh, Slopes S, Time Constants Tau
Vhm = -35*mV
Sm = -7.8*mV
Vhh = -55*mV
Sh = 7*mV
Vhn = -28*mV
Sn = -15*mV
VhmN = -30*mV
SmN = -5*mV
VhhN = -45*mV
ShN = 5*mV
VhmL = -40*mV
SmL = -7*mV
Vhmnap = -47.1*mV
Smnap = -3.1*mV
Vhhnap = -59*mV
Shnap = 8*mV
TaumN = 4 *ms
TauhN = 40*ms
TaumL = 40*ms
Taum = 0.000001*ms
Taumnap = Taum

# Reversal potentials in mV
ENa = 55 * mV
EK = -80 * mV
ECa = 80 * mV
Eleak = -60 *mV



def iext1(t):
    if t <= 500*ms:
        return 0*uA/cm**2
    elif t <= 5500*ms:
        return (t/ms-500) * 25/5000 *uA/cm**2
    elif t <= 15000*ms:
        return (25-(t/ms-5500) * 25/5000) *uA/cm**2
    else:
        return 0 *uA/cm**2
def iext(t):
    if t <=500*ms:
        return -20*uA/cm**2
    elif t <= 5500*ms:
        return (-20+(t/ms-500) * 25/5000) *uA/cm**2
    elif t <= 100000*ms:
        return (5-(t/ms-5500) * 25/5000) *uA/cm**2
    else:
        return 0


eqsMN=Equations('''
dh/dt= (hinf-h)/Tauh : 1
dn/dt= (ninf-n)/Taun : 1
dmnS/dt = (mnSinf-mnS)/TaumN : 1
dhnS/dt = (hnSinf-hnS)/TauhN : 1
dmnD/dt = (mnDinf-mnD)/TaumN : 1
dhnD/dt = (hnDinf-hnD)/TauhN : 1
dmnap/dt = (mnapinf-mnap)/Taumnap : 1
dhnap/dt = (hnapinf-hnap)/Tauhnap : 1
dm/dt = (minf-m)/Taum : 1
dml/dt = (mlinf-ml)/TaumL : 1
dCaS/dt = f*(-alpha*ICaS-kca*CaS) : mole*dmetre**-3
dCaD/dt = f*(-alpha*ICaD-kca*CaD) : mole*dmetre**-3
dVs/dt = 1/C*(iext(t)-INaS-IKS-ICaS-IleakS-IcouplingS) : mV
dVd/dt = 1/C*(-INapD-IKD-ICaD-IleakD-IcouplingD) : mV

Tauh = 30/(exp((Vs/mV+50)/15)+exp(-(Vs/mV+50)/16))*ms : ms
Taun = 7/(exp((Vs/mV+40)/40)+exp(-(Vs/mV+40)/50))*ms : ms
Tauhnap = 1200/(cosh(Vd/mV + 59)/16)*ms :ms
minf = 1/(1+exp((Vs-Vhm)/Sm)) : 1
hinf = 1/(1+exp((Vs-Vhh)/Sh)) : 1
ninf = 1/(1+exp((Vs-Vhn)/Sn)) : 1
mnapinf = 1/(1+exp((Vd-Vhmnap)/Smnap)) : 1
hnapinf = 1/(1+exp((Vd-Vhhnap)/Shnap)) : 1
mnSinf = 1/(1+exp((Vs-VhmN)/SmN)) : 1
hnSinf = 1/(1+exp((Vs-VhhN)/ShN)) : 1
mnDinf = 1/(1+exp((Vd-VhmN)/SmN)) : 1
hnDinf = 1/(1+exp((Vd-VhhN)/ShN)) : 1
mlinf = 1/(1+exp((Vd-VhmL)/SmL)) : 1

INapD = GNapD*mnap*hnap*(Vd-ENa) : mA*umetre**-2
INaS = GNa*m**3*h*(Vs-ENa) : mA*umetre**-2
IKS = (GK_dr*n**4 + GK_CaS*CaS/(CaS+Kd))*(Vs-EK) : mA*umetre**-2
ICaS = GCa_NS*mnS**2*hnS* (Vs-ECa) : mA*umetre**-2
IleakS = gleak*(Vs-Eleak) : mA*umetre**-2
IcouplingS = gc/p*(Vs-Vd) : mA*umetre**-2
IKD = GK_CaD*CaD/(CaD+Kd)*(Vd-EK) : mA*umetre**-2
ICaD = (GCa_ND*mnD**2*hnD+GCa_L*ml)*(Vd-ECa) : mA*umetre**-2
IleakD = gleak*(Vd-Eleak) : mA*umetre**-2
IcouplingD = gc/(1-p)*(Vd-Vs) : mA*umetre**-2
''')

P=NeuronGroup(1,model=eqsMN,
threshold=EmpiricalThreshold(threshold=-40*mV),
implicit=True)

ivs= -60
ivd= -60
P.Vs= ivs*mV
P.Vd= ivd*mV
#P.h = 1/(1+exp((ivs-Vhh)/Sh))
P.h = 0.9
#P.n = 1/(1+exp((ivs-Vhn)/Sn))
P.n = 0
#P.mnS = 1/(1+exp((ivs-VhmN)/SmN))
P.mnS = 0
#P.hnS = 1/(1+exp((ivs-VhhN)/ShN))
P.hnsS = 0.9
#P.mnD = 1/(1+exp((ivd-VhmN)/SmN))
P.mnD = 0
#P.hnD = 1/(1+exp((ivd-VhhN)/ShN))
P.hnD = 0.9
#P.ml = 1/(1+exp((ivd-VhmL)/SmL))
P.ml = 0
P.m = 0
P.CaS = 0
P.CaD = 0

tracevm=StateMonitor(P,'Vs',record=[0])
traceVd=StateMonitor(P,'Vd',record=[0])
spikes = SpikeMonitor(P)

figure(1)
ion()
run(1*ms)
subplot(211)
tracevm.plot(refresh=100*ms,showlast=15000*ms)
subplot(212)
traceVd.plot(refresh=100*ms,showlast=15000*ms)
run(14999*ms)
ioff()
show()